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Abstract 

MicroRNAs (miRNAs) are involved in plant growth, development, and response to biotic and abiotic stresses. Most of the 
miRNAs that have been identified in model plants are well characterized, but till now, no reports have previously been 
published concerning miRNAs in Paulownia australis. In order to investigate miRNA-guided transcript target regulation in P. 
australis, small RNA libraries from two P. australis (diploids, PA2; and autotetraploids, PA4) genotypes were subjected to 
Solexa sequencing. As a result, 10,691,271 (PA2) and 10,712,733 (PA4) clean reads were obtained, and 45 conserved miRNAs 
belonging to 15 families, and 31 potential novel miRNAs candidates were identified. Compared with their expression levels 
in the PA2 plants, 26 miRNAs were up-regulated and 15 miRNAs were down-regulated in the PA4 plants. The relative 
expressions of 12 miRNAs were validated by quantitative real-time polymerase chain reaction. Using the degradome 
approach, 53 transcript targets were identified and annotated based on Gene Ontology and Kyoto Encyclopedia of Genes 
and Genomes pathway analysis. These targets were associated with development, stimulus response, metabolism, signaling 
transduction and biological regulation. Among them, 1 1 targets, including TCP transcription factors, auxin response factors, 
squamosa promoter-binding-like proteins, scarecrow-like proteins, L-type lectin-domain containing receptor kinases and 
zinc finger CCCH domain-containing protein, cleaved by four known miRNA family and two potentially novel miRNAs were 
found to be involved in regulating plant development, biotic and abiotic stresses. The findings will be helpful to facilitate 
studies on the functions of miRNAs and their transcript targets in Paulownia. 
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Introduction 

MicroRNAs (miRNAs), a class of endogenous 21-24 nucleotide 
(nt), single-stranded, non-coding RNAs, modulate the gene 
expression at the transcription and post-transcriptional levels 
[1,2]. In plants, emerging data demonstrate that miRNAs play 
vital regulatory roles in a wide range of biological processes, 
including regulation of plant growth, development, signal trans- 
duction [3-6], and response to biotic and abiotic stresses via 
interactions with their specific target mRNAs [7—10]. Thereby, 
understanding of miRNAs functions in plant requires the 
recognition of their target genes. Initially, based on the either 
perfect or near-perfect sequence complementarity between 
miRNAs and their target mRNAs, the plant miRNA transcript 
targets have been identified via computational prediction [11,12]. 
However, computational target prediction method is very 
challenging because of the existence of a higher mismatch in 
miRNA-target mRNA pairing [13]. With the development of 
high-throughput sequencing technology, recently, a transcriptome 
wide experimental method combining high-throughput miRNA 



profiling with degradome sequencing analysis, has been success- 
fully developed to identify miRNA-directed mRNA cleavage at a 
large scale [14,15]. Using this strategy, so many plant species have 
been studied, for instance, Arabidopsis [15], Oiyza saliva [16], 
maize [17], Brassica napus [18], Medicago truncatula [19], 
Fragaria ananassa [20], grapevine [21], Taxus [22], Populus 
trichocarpa [23], Populus euphratica [24], Trifoliate orange [25], 
Paulownia tomentosa [26], and Paulownia fortunei [27]. These 
discoveries have triggered detailed biological studies on gene 
regulation with the result that the number of reported plant 
miRNAs and their transcript targets has been increasing rapidly. 

Paulownia is a genus of deciduous tree species in the family 
Paulowniaceae, related to and sometimes included in the family 
Scrophulariaceae. It is a very adaptable, extremely fast-growing 
woody plant, and in managed plantations, it can be harvested for 
saw timber in as little as 5 years. Paulownia can not only provide 
wood for a variety of products including timber, fuelwood, herbal 
medicines, boxes, clogs, musical instruments, and surfboards, but 
can also benefit the environment and ecology; for example, it can 
increase crop production when used for intercropping on 
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farmland, and prevent soil erosion [28,29] . Because of its specific 
characteristics and economic value, in its native China, Paulownia 
is popular for reforestation, intercropping on farmland, and 
roadside and garden planting. To enlarge the germplasm and 
breed new varieties, in recent years, several autotetraploid 
paulownia seedlings were induced successfully using colchicines 
[30-33]. The autotetraploid cultivars of paulownia contains two 
sets of the same chromosomes, and compared with their diploid 
progenitors, the autotetraploids have better timber quality and 
improved stress resistance [34,35], which makes understanding the 
molecular mechanisms underlying the differences of diploids and 
autotetraploids imperative. To our knowledge, the different 
paulownia species have genetic diversity and differentiation, 
although miRNAs in the P. tomentosa [26] and P. fortunei 
species have been identified [27], there have been no detailed 
reports of miRNAs in P. australis. In the current study, we used 
Solexa sequencing and degradome technology to analyze four 
sequencing libraries constructed from the seedlings of the two P. 
australis genotypes to identify conserved and novel miRNAs and 
to investigate the potential roles of their transcript targets. 

Materials and Methods 

Plant culture and treatment 

Plantlets of two P. australis genotypes (diploids, PA2; and 
autotetraploids, PA4) were cultured in vitro on medium as 
described previously [30,36]. All the samples were grown at 
25°C under a 16/8 h (light/dark) photoperiod for 30 days, and at 
least three parallel samples were prepared for each genotype. After 
culturing, samples were collected from nine different seedlings of 
each genotype and subsequently mixed, frozen immediately in 
liquid nitrogen, and stored at — 86°C until further use. 

Construction and sequencing of the small RNA library 

Total RNA was isolated from the leaves of the PA2 and PA4 
samples using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) 
following the manufacturer's instructions. Two small RNA (sRNA) 
libraries were constructed, one from the PA2 sample and the other 
from the PA4 sample, using methods described elsewhere [24,37]. 
Briefly, bands in the 18-30 nt size range were purified by 
electrophoretic separation on a 15% denaturing polyacrylamide 
gel and then successively ligated with 5' and 3' adapters. After 
reverse transcription and PCR, the amplified products were 
subjected to Solexa sequencing (Illumina, San Diego, CA, USA) at 
the Beijing Genomics Institute, Shenzhen, China. The data used 
in this publication have been deposited in the NIH Short Read 
Archive database (http://www.ncbi.nlm.nih.gov/sra) with acces- 
sion number SRP041440 and SRP041442 (Alias: PRJNA245384 
and PRJNA245383). 

Identification of conserved and novel miRNAs in the P. 
australis libraries 

Low-quality reads, poly(A) reads, oversized insertions, reads 
shorter than 18 nt, and adaptor contaminated reads were filtered 
out from the raw reads to yield the clean reads. The length 
distribution of the 18-30 nt clean reads was analyzed and the 
reads were mapped to the sequences in the P. australis UniGene 
database (http://www.ncbi.nlm.nih.gov/sra; Accession Number 
SRP032321) using SOAP [38]. Perfecdy matched reads were 
analyzed by running Blastall (http://www.ncbi.nlm.nih.gov/staff/ 
tao/URLAPI/blastall/) against the Rfam (http://www.sanger.ac. 
uk/software/Rfam) and GenBank databases (http://www.ncbi. 
nlm.nih.gov/) to discard tRNAs, rRNAs, snRNAs, snoRNAs, and 
other ncoRNAs. The remaining reads were annotated by aligning 



them against the sequences in the miRBase database 19.0 (http:// 
www.mirbase.org/) to identify known miRNAs. When a read in 
the P. australis libraries shared homology with fewer than two 
mismatches with a miRNA sequence in miRBase, it was 
considered to be an evolutionarily conserved miRNA. 

The MIREAP software (https://sourceforge.net/projects/ 
mireap/) was used to detect novel miRNAs in the two libraries 
by predicting the stem-loop structures and estimating the minimal 
folding free energy of the unannotated sRNAs (no matches in 
miRBase) that mapped to the P. australis UniGene sequences. 
The stem-loop structures of the candidate precursor miRNAs (pre- 
miRNAs) were predicted using Mfold (http://mfold.rna.albany. 
edu/?q0mfold) [39]. The strict criteria used to annotate the 
candidate miRNAs were as described by Meyers el al. [40]. 

Differential expression analysis of miRNAs in the two P. 
australis libraries 

To identify the differentially expressed miRNAs between the 
PA2 and PA4 libraries, the abundance of the miRNAs in the two 
libraries were normalized to one million, regardless of the total 
number of miRNAs in each sample. Then, statistical analysis was 
performed based on a Poisson distribution. Finally, when the fold- 
change was greater than 1.0 or less than —1.0 and the P-values 
were less than 0.05, the miRNA was considered to be significandy 
differentially expressed. The fold-changes and P- values [41,42] 
were calculated from the normalized expression using the 
following formulas: 

Normalized expression = (actual miRNA reads/ total count of 
clean reads) x 1,000,000 

Fold-change = log2 (PA4 normalized reads/PA2 normalized 
reads) 

P -values were calculated as, 
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Degradome library construction, data analysis, and target 
identification 

Two degradome libraries were constructed from the leaves of 
the PA2 and PA4 plants based on a protocol published previously 
[14,15]. Briefly, using T4 RNA ligase (Takara, Dalian, China), 
poly(A) enriched RNA was ligated to a 5 'RNA adapter containing 
a Mmel recognition site. Reverse transcription was performed to 
generate first-strand cDNA, followed by PCR enrichment, and 
digestion using the Mmel restriction enzyme (NEB, Ipswich, MA, 
USA). A double-stranded DNA adapter was then ligated to the 
digested products using T4 DNA ligase (NEB, Ipswich, MA, USA) 
and gel purified for PCR amplification. The resulting PCR 
products were sequenced using an Illumina HiSeqTM 2000 
system. Low-quality sequences and adapters were removed from 
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the raw reads of the degradome libraries, then the remaining reads 
were aligned to the P. australis UniGene database (http:/ /www. 
ncbi.nlm.nih.gov/sra; Accession Number SRP032321) using 
SOAP software (http://soap.genomics. org.cn/) to define the 
coverage rate. The PAIRFINDER software (version 2.0) [22] was 
used to identify the miRNA-mediated cleaved fragments. Align- 
ments with scores not exceeding four and no mismatches at the 
cleavage site (between the 10th and 11th nucleotides) were 
considered to be the potential miRNA targets. Furthermore, t- 
plots were built according to the distribution of signatures (and 
abundances) along these transcripts to help analyze the miRNA 
targets and RNA degradation patterns. Based on the locations of 
the target genes in the P. australis UniGene sequences, putative 
target genes were selected manually and subsequently mapped to 
the previously identified genes that were annotated according to 
the annotations of their homologous sequences in the GenBank Nr 
and Nt databases, and in the Swiss-Prot database using BLASTX 
searches with an E-value cutoff of 10 5 . 

Verification of miRNAs and their targets by quantitative 
real-time polymerase chain reaction 

Quantitative real-time polymerase chain reaction (qRT-PCR) 
was used to validate and measure the expressions of selected 
miRNAs and their targets obtained from the Solexa sequencing 
and degradome analysis. The qRT-PCR was performed according 
to the protocol described previously [43] . For the qRT-PCR, total 
RNA were isolated from the leaves of two individual PA2 and PA4 
plants at two developmental stages (30-day-old in vitro plandets, 
and two-year-old saplings from field). Three biological replicates of 
each stage were used. The qRT-PCRs were performed in 
triplicate on a CFX96TM Real-Time PCR System (Bio-Rad, 
Hercules, CA, USA) using a Superscript III platinum SYBR 
Green one-step qRT-PCR kit (Invitrogen, Carlsbad, CA, USA) 
according to the manufacturer's instructions. The PCR conditions 
were 50°C for 3 min, 95°C for 5 min, then 40 cycles of 95°C for 
15 s, 55°C for 30 s, and 40°C for 10 min. Specific stem-loop 
primers and other primers for the miRNAs were designed based 
on the mature miRNA sequence. The U6 was used as an 
endogenous control. All the primers used for the qRT-PCR are 
listed in Table S 1 . The primers for the target genes were designed 
using Beacon Designer, version 7.7 (Premier Biosoft International, 
Ltd., Palo Alto, CA, USA), and the 18S rRNA of Paulownia was 
chosen as an endogenous reference gene for normalization. The 
relative expression level of a miRNA or target gene was calculated 
according to the method described previously [44] . 

Results 

Analysis of sRNAs in the two P. australis libraries 

A total of 13,895,340 and 13,537,466 raw reads were obtained 
in the PA2 and PA4 sRNA libraries, respectively. After removing 
low-quality reads, poly(A) reads, oversized insertions, reads shorter 
than 18-nt, and adaptor contaminated reads, we obtained 
10,691,271 (PA2) and 10,712,733 (PA4) clean reads that 
represented 2,006,153 (PA2) and 2,418,971 (PA4) unique reads 
with lengths ranging from 18 to 30-nt. The majority of clean reads, 
63.24% (PA2) and 68.79% (PA4), ranged from 20 to 24-nt in 
length, among which the 24-nt long reads were the most 
abundant, followed by the 21-nt long reads (Figure SI). The 
abundances of the 24-nt and 21-nt long sRNAs in the PA4 library 
were 6.48% and 0.04% more, respectively, than their abundances 
in the PA2 library. This atypical situation has also been reported in 
other hardwood species, such as Liriodendron chinense [45]. 



The clean reads were used to query the P. australis UniGene 
database (http://www.ncbi.nlm.nih.gov/sra), the non-coding 
RNAs sequences deposited in GenBank (http://www.ncbi.nih. 
gov/Genbank), the Rfam database (http://rfam.sanger.ac.uk/), 
and miRBase 19.0 (http://microrna.sanger.ac.uk/sequences). 
These searches allowed us to assign annotations to each sRNA 
sequence, including rRNA, tRNA, snRNA, snoRNA, and known 
miRNA (Table S2). As a result, 848,652 (PA2) and 931,989 (PA4) 
known miRNA reads were detected. Finally, 1,809,478 (PA2) and 
2,308,527 (PA4) unique unannotated sRNAs remained for 
predicting potentially novel miRNAs. 

Identification of conserved miRNAs in P. australis 

To detect conserved miRNAs in the two libraries, the unique 
reads were compared with the known miRNAs in miRBase 19.0, 
allowing two mismatches. In the two libraries, a total of 45 
conserved miRNAs belonging to 15 miRNA families were found 
to share high identity with known plant miRNAs (Table 1 and 
Table S3). Among these 45 miRNAs, 40 and 43 were identified in 
PA2 and PA4, respectively; 38 of them were expressed in both 
libraries, two were found only in the PA2 library, and five 
occurred only in the PA4 library. Most of the conserved miRNAs 
(91.67%) were 21-nt long and the remaining miRNAs were 20 or 
23-nt long (Table 1 and Table S3). This result is consistent with 
the current understanding that canonical miRNAs are 2 1 -nt long, 
while canonical small interfering RNAs (siRNAs) are 24-nt long 
[46]. Among the 15 miRNA families, the miR166 family was 
dominant in both libraries, accounting for 89.44% (PA2) and 
86.35% (PA4) of all conserved miRNA reads, followed by the 
miR159 family (Table 1). Several miRNA families, such as 
miR156, miR396, miR397, miR398, and miR482, had moderate 
expression levels in the two libraries, while other miRNA families 
showed very low levels of expression, with fewer than 100 reads, in 
both libraries. Moreover, different members in the same miRNA 
family displayed significantly different expression levels. For 
instance, members of the miR166 family varied in abundance 
from 50,932 to 563,946 reads in the two libraries (Table 1). 

Identification of novel miRNAs in the P. australis sRNA 
libraries 

The characteristic stem-loop structure of pre-miRNA was 
employed to predict novel miRNAs [40]. We found that 31 of 
the potential pre-miRNAs met this requirement; 15 pre-miRNAs 
exhibited the 3p:5p miRNA pair, providing more evidence that 
they were novel miRNAs, while the other 16 miRNAs were also 
considered to be potentially novel miRNAs (Table 2 and Table 
S4). Among the 31 candidate novel miRNAs, 10 miRNAs were 
common between the two libraries, while 10 and 1 1 were specific 
to PA2 and PA4, respectively. The length of the mature miRNAs 
varied from 20 to 23 nt, with the majority being 21-nt long. The 
mature miRNA sequences were localized inside the stem-loop 
structures, with almost half in either the 3p or 5p arms. We 
observed that the average pre-miRNA length was 151 nt, and 
these precursors had minimal folding free energies that ranged 
from —142.5 kcal mol -1 to —30.8 kcal mol 1 with an average of 
— 56.3 kcal mol . In this study, compared with the abundance of 
the conserved miRNAs, the majority of novel miRNAs had 
relatively low expressions, and only five (PA2) and six (PA4) of the 
novel miRNA candidates had more than 1,000 reads (Table 2 and 
Table S4). 
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Expression patterns of conserved and novel miRNAs 

The differential expression analysis was performed based on the 
normalized read counts for each identified miRNA. A total of 41 
miRNAs showed statistically significant changes (fold change £ 1 .0 
or < — 1 .0, and P-values ^0.05) in their relative abundance 
between the two libraries. Among them, 26 miRNAs were up- 
regulated (13 conserved miRNAs and 13 novel miRNAs; fold- 
change S1.0) and 15 were down-regulated (three conserved 
miRNAs and 12 novel miRNAs; fold-change <-1.0) in the PA4 
library compared with in the PA2 library. The expression levels of 
some miRNAs changed significantly. For example, the expression 
levels of pas-miR169b-3p, pas-miR169c-3p, pas-miR396c-3p, pas- 
miR396d-3p, pas-miR171a, pas-miR171b, and pas-miR171c 
increased or decreased by about 5-fold in the PA4 library 
(Table 1). We also detected 21 novel miRNAs in the PA4 library 
that had 5-fold greater or lesser expression levels than in the PA2 
library (Table 2). 
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Identification of miRNAs transcript targets in P. australis 
by using degradome analysis 

To better understand the physiological functions and biology 
processes that these miRNAs may be involved in during the 
development of P. australis, transcript target identification was 
performed based on the degradome approach. Using PAIRFIN- 
DER, a total of 53 targets were predicted to be cleaved by 1 1 of 
the conserved miRNA families and three novel miRNA candidates 
(Table S5). The target transcripts were pooled and grouped into 
three categories (I — III) according to their relative abundances 
[14,15] (Figure 1). Among these identified targets, 32 (37 cleavage 
sites) belonged to category I, and 20 (34 cleavage sites) and six (six 
cleavage sites) targets belonged to categories II and III, respec- 
tively (Table S5). We then performed a Gene Ontology (GO) 
analysis to assign functional annotations to the predicted target 
genes, as described previously [47]. More than 80% of the target 
genes were annotated as being involved in regulation of biological 
processes, and the GO annotations for the predicted target genes 
are shown in Figure 2. The Kyoto Encyclopedia of Genes and 
Genomes (KEGG) Pathway database was used to further classify 
the miRNA target genes. The KEGG Pathway annotation showed 
that the target genes were involved in metabolic pathways, plant- 
pathogen interaction, pyruvate metabolism, carbon fixation in 
photosynthetic organisms, purine metabolism, RNA polymerase, 
pyrimidine metabolism, nitrogen metabolism, plant hormone 
signal transduction, cellular metabolism, and disease. (Table S5). 
The prediction and annotation of the miRNA target genes may 
provide some new insights into how P. australis miRNAs regulate 
gene expression in this plant. 
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Validation of the candidate miRNAs and their targets by 
qRT-PCR 

To verify the existence and expression levels of the miRNAs 
determined from the high-throughput Solexa sequencing, 12 
miRNAs with different expression levels were randomly selected 
for qRT-PCR analysis. As shown in Figure S2, the expression 
patterns of the miRNAs obtained by qRT-PCR showed similar 
trends to their expression patterns in the two libraries as 
determined by Solexa sequencing. Compared with their expres- 
sions in PA2, the expression levels of pas-miR156c, pas-miR398a- 
3p, pas-miR408a-5p, and pas-mir22-3p were up-regulated in PA4 
in the 30-day plantlets and in the 2-year saplings, pas-miR3 1 9a-3p 
were on the contrary. Further, the expressions of pas-miR 160a, 
pas-miR167, pas-miR171a, pas-miR397a, pas-mirl, and pas- 
mirl4 in PA4 were up-regulated in the 30-day plantlets, and 
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Table 2. Novel miRNAs identified from P. australis. 



miRNA 


expression 




Fold-change 


P-value 


MiRNA* 


expression 




PA2 


PA4 


(log2PA4/PA2) 




PA2 


PA4 


pas-mirl 


116 


422 


1.86 


7.02E-42 


106 


74 


pas-mir2 


43 


0 


-8.65 


1.09E-13 


0 


0 


pas-mir3 


8 


3 


-1.42 


1.45E-01 


6 


5 


pas-mir4a 


6 


0 


-5.81 


1.55E-02 


0 


0 


pas-mir4b 


6 


0 


-5.81 


1.55E-02 


0 


0 


pas-mir5a 


7116 


6078 


-0.23 


5.38E-20 


208 


223 


pas-mir5b 


7116 


6078 


-0.23 


5.38E-20 


208 


223 


pas-mir6a-3p 


2125 


4098 


0.94 


1.51E-139 


67 


71 


pas-mir6b-3p 


2125 


4098 


0.94 


1.51E-139 


67 


71 


pas-mir6c-3p 


2125 


4098 


0.94 


1.51E-139 


67 


71 


pas-mir7-3p 


696 


0 


-12.67 


1.51E-210 


6 


0 


pas-mir8a-3p 


29 


0 


-8.08 


4.57E-01 


0 


0 


pas-mir8b-3p 


29 


11 


-1.4 


5.48E-02 


0 


0 


pas-mir9 


8 


0 


-6.23 


3.87E-03 


0 


0 


pas-mirl 0-3p 


10 


0 


-6.55 


9.66 E-04 


0 


0 


pas-mirl 1-3p 


17 


0 


-7.31 


7.49 E-06 


0 


0 


pas-mirl 2 


5 


0 


-5.55 


3.1 1 E-02 


0 


0 


pas-mirl 3-3p 


934 


1506 


0.69 


5.10E-31 


1128 


1403 


pas-mirl 4 


10 


28 


1.48 


3.45 E-03 


2 


1 


pas-mirl 5a 


345 


0 


-11.66 


9.86E-105 


107 


0 


pas-mirl 6a-3p 


0 


19 


7.47 


1.95 E-06 


0 


0 


pas-mirl 6b-3p 


0 


19 


7.47 


1.95 E-06 


0 


0 


pas-mirl 6c-3p 


0 


19 


7.47 


1.95 E-06 


0 


0 


pas-mirl 7-3p 


0 


8 


6.22 


3.94E-03 


0 


7 


pas-mirl 8-3p 


0 


11 


6.68 


4.94E-04 


0 


0 


pas-mirl 9-3p 


0 


126 


10.2 


1 .34E-38 


0 


0 


pas-mir20-3p 


0 


54 


8.98 


5.87E-17 


0 


2 


pas-mir21a 


0 


12 


6.81 


2.47E-04 


0 


9 


pas-mir21b 


0 


12 


6.81 


2.47E-04 


0 


9 


pas-mir22-3p 


0 


11 


6.68 


4.94E-04 


0 


0 


pas-mir23-3p 


0 


28 


8.03 


3.84E-09 


0 


0 
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down-regulated in the PA4 in the 2-year saplings, while the 
expressions of pas-mir3 were the opposite to this. Thus, with 
the development of the plants, the expression levels of some of the 
miRNAs showed different trends. Six miRNAs (pas-miR 1 60a, pas- 
miR167, pas-miR319a-3p, pas-miR398a,-3p pas-miR408a-3p, 
pas-mirl) had the same expression trend between PA2 and PA4 
in the two stages. Furthermore, to confirm the reliability of 
degradome sequencing technology and the potential correlation 
between miRNAs and their transcript targets, 12 genes from two 
P. australis genotypes were also selected for qRT-PCR assays. 
The results showed that except the targets of auxin response factor 
(ARF) ARF8 (CL421 l.Contig3) and scarecrow-like protein (SCL) 
SCL15 (CL10503.Contigl), the expression levels of the rest targets 
were inversely correlated with these of the corresponding miRNAs 
(Figure 3). During the different developmental stages, pas- 
miR160a, pas-miR167, pas-miR171 and pas-mirl at a relatively 
higher level in the PA4 than in the PA2 at the 30-day plantlets 
stage, and lower level at the 2-year saplings stage, while its 



transcript targets, CL3173.Contig7, CL1 1603.Contigl, 
CL6407.Contig9, CL1 1078.Contig2, CL1 1078.Contig3, and Uni- 
gene9061 expressed in the reverse way as expected, and these 
coding for proteins are members of the auxin response factors 
ARF 10, ARF 18 and ARF6, the scarecrow-like proteins SCL6 and 
SCL22, and serine/threonine protein kinase, respectively. The 
expression levels of pas-miR319a-3p in PA4 were significandy 
lower than in the PA2 at two treatment stages, while the reverse 
was true for the transcription factor TCP4 (CL9103.Contig3) 
(Figure 3). Moreover, a reverse trend was noted between pas- 
miR156c and its target genes coding for squamosa promoter- 
binding-like protein (SPL) SPL6 and SPL12 (CL1 1428.Contig2 
and CL5129.Contig2), and between pas-mir22-3p and its target 
gene coding for zinc finger CCCH domain-containing protein 53 
(CL1 197.Gontig2) in the PA4 as compared to the PA2 (Figure 3). 
These results indicated that the miRNA and transcript target 
expression patterns were very complex and varied during the 
growth and development of P. australis. Furthermore, the possible 
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Figure 1. Target plots (t-plots) of miRNA targets in different categories confirmed by degradome sequencing. (A) T-plot (top) and 
miRNA: mRNA alignments (bottom) for two category I targets, Unigenel 1776_AII and CL1057.Contig2_AII transcripts. The arrow indicates signatures 
consistent with miRNA-directed cleavage. The solid lines and dot in miRNA: mRNA alignments indicate matched RNA base pairs and GU mismatch, 
respectively, and the red letter indicates the cleavage site. (B) CL5129.Contig3_AII and CL6407.Contig9_AII, a category II target for pas-miR156d and 
pas-miR167. (C) CL1 1 907.Contig 1 _AII and Unigene1681_AII, a category III target for pas-miR168b and pas-miR169e. 
doi:10.1371/journal.pone.0106736.g001 



roles of these miRNAs in the genome duplication changes from 
diploid to autotetraploid in these plants were revealed. 

Discussion 

Variations in plant morphology and physiology resulting from 
genome duplication have occurred in many plants, such as 
Triticum, Gossypium, Spartina, Tragopogon, Brassica and Sola- 
rium [48,49]. Generally speaking, genome duplication has led to 
the production of fast-growing, high-quality plants [50,51]. The 
autotetraploid P. australis, which contains two sets of the same 
chromosomes, displayed apparent alterations in morphology, 
growth development, physiology, and gene expression when 
compared with their diploid counterparts [34,35]. However, the 
mechanisms for these changes are poorly understood. MiRNAs 
are a class of endogenous small RNAs that have been involved in 
many processes, including growth, development, and resistance to 
stress and disease, by their ability to regulate gene expression in 
plants [52-54]. To understand the functions of miRNAs in diploid 
and autotetraploid P. australis, in this present study, we used 
Solexa sequencing and degradome approaches to construct two 
sRNA libraries and two degradome libraries from the PA2 and 
PA4 plants to identify conserved and novel miRNAs and their 
transcript targets. A total of 45 conserved miRNAs belonging to 15 
miRNA families and 3 1 potential novel miRNA candidates along 
with 53 transcript targets were identified across the PA2 and PA4 
libraries. Most of the identified conserved miRNA families are also 
conserved in other plant species, including Populus tomentosa 
[53], Populus euphratica [55], Oryza saliva [3], and Arabidopsis 
thaliana [56]. In the PA2 and PA4 libraries, the expression 



patterns varied dramatically among the different miRNA families, 
and different members in the same miRNA family also displayed 
significandy different expression levels. For instance, the read 
number varied from 26 (miR171 family) to 1,422,532 (miR166 
family) (Table 1), and members of the miR166 family varied in 
abundance from 50,932 to 563,946 reads (Table 1). Moreover, the 
majority of novel miRNAs had relatively low expressions, and only 
five (PA2) and six (PA4) of the novel miRNA candidates had more 
than 1,000 reads (Table 2). Our results are in accordance with 
previous reports that novel miRNAs were often represented in 
relatively lower levels than conserved miRNAs [57,58]. Further- 
more, it is possible that the low-expression novel miRNAs may 
play particular functions in specific tissues, during developmental 
stages, or under various growth conditions. Whether these low- 
expression miRNAs are expressed at higher levels in other tissues 
or developmental stages, or are regulated by environmental 
stresses, remain to be investigated. 

Previous studies have shown that many of the genes appeared to 
be methylated in tetraploid Paulownia plants specifically after 
genome duplication by the methylation-sensitive amplified poly- 
morphism analysis [59]. DNA methylation has been reported to be 
involved in inducing gene silencing, which can restart or change 
the genes expression levels [60]. Salmon et al. [61] found that 
significant changes in DNA methylation patterns could explain the 
morphological plasticity and larger ecological amplitude of 
Spartina allopolyploids. Indeed, we found that the expression 
levels of many of the differentially expressed miRNAs in the PA4 
library were not increased by more than two-fold compared with 
their expressions in the PA2 library. However, the expressions of 




Figure 2. Gene Ontology analysis of miRNA targets in P. australis. 
doi:1 0.1 371 /journal.pone.01 06736.g002 
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Figure 3. Relative expression levels of the target genes in P. australis. PA2-1, 30-day-old diploid in vitro plantlets; PA4-1, 30-day-old 
autotetraploid in vitro plantlets; PA2-2, two-year-old diploid saplings; PA4-2, two-year-old autotetraploid saplings. CL1 1428.Contig2 (SPL6) and 
CL5129.Contig2 (SPL12) targeted by pas-miR156c, CL3173.Contig7 (ARF10) and CL1 1603.Contig1 (ARF18) targeted by pas-miR160a, CL421 1.Contig3 
(ARF8) and CL6407.Contig9 (ARF6) targeted by pas-miR167, CL1 1078.Contig2 (SCL6), CL10503.Contig1 (SCL15) and CL1 1078.Contig3 (SCL22) targeted 
by pas-miR1 71a, CL9103.Contig3 (Transcription factor TCP4) targeted by pas-miR319a-3p; Unigene9061 (Serine/threonine protein kinase) targeted by 
pas-mir1 ; CL1 1 97.Contig2 (Zinc finger CCCH domain-containing protein 53) targeted by pas-mir22-3p. Three independent biological replicates were 
performed. Values are means ± SD (n = 3).The expression levels of targets were normalized to 18SrRNA. The normalized miRNA levels in the PA2-1 
were arbitrarily set to 1 . *: Statistically significant differences between PA2 and PA4 under the same developmental stages (p-value was less than 0.05) 
doi:1 0.1 371 /journal.pone.01 06736.g003 



about half of these miRNAs were significantly different in the two 
libraries. Some of the miRNAs were expressed at similar level in 
the PA4 and PA2 libraries (Table 2). These findings are similar to 
those reported previously in P. tomentosa, P. fortunei and 
tetraploid Arabidopsis thaliana lines [26,27,62], suggesting that 
the genome merger in the PA4 plants lead to nonadditive 
expression of the miRNA primary transcripts and miRNA target 
genes. Furthermore, the expression patterns of the differentially 
expressed miRNAs and transcript targets at different development 
stages were validated by qRT-PCR. The result showed that the 
differentially expressed miRNAs caused different expression levels 
in their transcript targets. Interestingly, we also found the 
expression levels of two target genes (CL421 l.Contig3 and 
CL10503.Contigl) were inconsistent with those of their corre- 
sponding miRNAs (pas-miR167 and pas-miR171a). The similar 
phenomena were also observed in P. tomentosa, Phalaenopsis 
Aphrodite and cotton [26,63,64], indicating that the other 
mechanisms of regulating expression of the target genes exist. 
Above all, these results imply that the miRNAs with significantly 
varied expressions in the PA2 and PA4 P. australis plants are 
probably involved in the epigenetic changes of PA4 plants; 
however, the relation between the miRNA expression patterns and 
genome duplication may be more complex than we first thought. 

To understand the biological functions of miRNAs, it is 
necessary to identify their transcript targets. In the present study, 
to avoid false-positive predictions of miRNA transcript targets in 
P. australis, we identified the 53 transcript targets for 1 1 miRNA 
families and three novel miRNA candidates by degradome 
sequencing, which opens up a new avenue for high-throughput 
validation of splicing targets [14,15]. The target genes predicted 
for the conserved miRNAs in P. australis were similar or 
functionally related to validated plant miRNA targets, which were 
annotated as being involved in diverse physiological processes. For 
instance, pas-miR156 targeted the SPL protein family, which can 
affect diverse developmental processes such as leaf development, 
shoot maturation, phase change, and flowering in plants [65]; and 
pas-miR167 targeted ARF6 and ARF 8, which belong to a class of 



transcription factors known to control multiple processes in plants, 
including the regulation of gynoecium and stamen maturation, 
and seed dispersal [66,67]. Thus, our results support the idea that 
conserved miRNAs take part in essential physiological processes in 
plants. 

The analysis of the target genes identified for the differentially 
expressed miRNAs revealed that some of the target genes may 
play important roles in plant morphology and physiology. We 
found that the expression level of pas-miR319a-3p decreased by 
about 3-fold in PA4 compared with its level in PA2, and pas- 
miR319a-3p was predicted to target TCP transcription factors, 
which are plant-specific transcription factors that have been shown 
to participate in specifying plant morphological traits, such as 
organ border delimitation [68], cell division and proliferation [69], 
flower and leaf shape, and shoot outgrowth [70]. In Arabidopsis , 
the TCP transcription factors have been related to control of the 
morphology of shoot lateral organs and formation of the shoot 
meristem-dependent regulation of the expression of boundary- 
specific genes [71]. The putative transcription factor SCL is a 
number of GRAS protein family that are involved in several 
aspects of plant growth and development, including control of 
asymmetric cell division, maintenance of stem cell status, and 
induction of the regeneration of the root tip after laser ablation, 
and SCL expression has been associated with auxin distribution in 
root apical meristem [72-74]. In this study, three SCL genes were 
predicted to be targeted by the miR171 family, which was up- 
regulated in the PA4 plants. This result suggested that SCL may 
be involved in the formation of adventitious root and the other 
aerial organs of Paulownia. L-type lectin-domain containing 
receptor kinases (lecRKs) and zinc finger CCCH domain- 
containing protein were predicted to be targeted by pas-miRl 
and pas-miR22-3p, respectively, suggesting that, besides their 
possible involvement in plant development, some of the miRNA 
target genes identified in this study could play fundamental roles in 
biotic and abiotic stresses. A previous study showed that LecRKs 
were most likely plasma membrane proteins, and were probably 
involved in mediating protein-protein interactions with a wide 
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range of functions such as recognition of oligosaccharide or 
lipochitooligosaccharide signals, linking ABA-signaling, and re- 
sponse to salt, drought, cold stress, wounding, and disease in plants 
[75,76]. Zinc finger GGCH domain-containing protein is a kind of 
RNA-binding protein and many studies have shown that it may be 
regulated by abiotic or biotic stresses, and could have regulatory 
functions in mRNA processing [77-80], thus supporting the 
possible roles for pas-miRl and pas-miR22-3p in the adaptive 
response of PA4 to abiotic stress. 

Based on the GO analysis of the targets of the identified 
miRNAs in the two P. australis genotypes plants, the target genes 
were separated into the biological process, cellular component, 
and molecular function. Some of the genes were annotated as 
being involved in biological regulation, cellular process, develop- 
mental process, response to stimulus, metabolic process, cell, cell 
part, organelle components, binding, and catalytic activity 
(Figure 2). Thus, our results suggest that the transcript targets 
might be closely related to the observed differences of phenotype 
(environmental adaptations) and resistance to biotic and abiotic 
stresses between the PA4 and PA2 plants. The functional role of 
the differentially expressed miRNAs will be the focus of future 
investigations. No target genes were predicted for many of the 
miRNAs identified in the two P. australis sRNA libraries, and a 
few of the predicted target genes were annotated as being of 
unknown function and hypothetical genes. Careful analysis of 
these potential targets will contribute further to our understanding 
of the role of miRNAs in P. australis. 

In summary, miRNA target genes were identified using a 
degradome approach that included functional annotation and 
pathway analyses based on the GO and KEGG databases. Some 
of the transcript targets regulated by the differentially expressed 
miRNAs were related to the physiology and environmental 
adaptations. Our results suggest that the significantly varied 
expression miRNAs in the two P. australis genotypes are probably 
involved in the epigenetic changes of PA4 plants, the correlation 
between the miRNA expression pattern and genome duplication 
may be more complex than we first thought. Taken together, this 
study provides beneficial information for elucidating the miRNA- 
mediated regulation of transcript targets in P. australis and related 
species. 
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